In the present specification, reference is made to the following publications illustrating prior art of conventional image reconstruction and imaging techniques.    [1] David Donoho. For most large underdetermined systems of linear equations, the minimal l1 norm near-solution approximates the sparsest near-solution. Communications on Pure and Applied Mathematics, 59(7):907-934, 2006;    [2] Jarvis Haupt and Robert Nowak. Signal reconstruction from noisy random projections. IEEE Trans. on Information Theory, 52(9):4036-4048, 2006;    [3] Maxim Raginsky, Rebecca M. Willett, Zachary T. Harmany, and Roummel F. Marcia. Compressed sensing performance bounds under poisson noise. IEEE Trans. on Signal Processing, 58(8):3990-4002, 2010;    [4] Emmanuel Candès and Justin Romberg. Sparsity and incoherence in compressive sampling. Inverse Problems, 23(3):969-985, 2007;    [5] C. C. Craig. On the Tchebychef inequality of Bernstein. Ann. Math. Statist., 4(2):94-102, 1933;    [6] J. M. Noras. Some formulas for moments of the Poisson distribution. Phys. Rev. B, 22(12):6474-6475, 1980;    [7] Emmanuel J. Candès et al. “An Introduction To Compressive Sampling” in “IEEE Signal Processing Magazine” March 2008, p. 21-30; and    [8] Rebecca M. Willett et al. “Performance bounds on compressed sensing with Poisson noise” in “IEEE International Symposium on Information Theory ISIT” 2009, Jun. 28, 2009-Jul. 3, 2009, p. 174-178; and    [9] PCT/EP2010/000526 (not published at the priority date of the present specification).
The reconstruction of a three-dimensional CT image is often accomplished by the well-known filtered back projection (FBP) algorithm, which is well suited for conventional CT imaging. However, this algorithm has some drawbacks. First, it requires an ideal measurement, i.e. the angles under which the patient is viewed must be many, and they should be equidistant. Second, the FPB does not work easily with noise. The latter means that the image acquisition must be done in such a way as to minimize noise, which requires a high dose. An alternative to FBP is to reconstruct the image iteratively using standard techniques such as Maximum Likelihood Expectation Maximization (MLEM). This method has the advantage that the imaging setup is in principle relatively arbitrary: for example, angles do not have to be equidistant (although they often are). Furthermore, the Poissonian nature of the noisy measurements is built into the method. This type of algorithm is becoming more and more feasible with the amount of computer power currently available. Nevertheless, it is still necessary to record a lot of information in order to obtain a good reconstruction in the end.
In contrast to this, the concept of compressive sensing (CS) promises, in principle, that a similarly good reconstruction can be obtained with significantly less recorded information. It relies on the idea that most, if not all, real world objects are “compressible” if they are represented in a suitable basis system. Compressible means that the coefficients of the representation, when ordered by size, fall off fast with some power of their index. The information that an object is compressible can be used in addition to the relatively few recorded data points to obtain a reliable reconstruction. With more details, compressed sensing (or: compressive sensing, compressive sampling and sparse sampling), is a technique for acquiring and reconstructing a signal utilizing the prior knowledge that it is sparse or compressible. An introduction to CS has been presented by Emmanuel J. Candès et al. [7]. The CS theory shows that sparse signals, which contain much less information than could maximally be encoded with the same number of data entries, can be reconstructed exactly from very few measurements in noise free conditions.
According to Jarvis Haupt et al. [2], a set of data f*j of size v (j=1, . . . , v) which is sparse (or rather “compressible”) can be accurately reconstructed from a small number k of random projections even if the projections are contaminated by noise of constant variance, e.g. Gaussian noise. Specifically, yi=Σjφijf*j+ξi with i=1, . . . , k are the noisy projections of f*j, taken with the projection matrix φij which consists of random entries all drawn from the same probability distribution with zero mean and variance 1/v (such that the transformed values yi have the same order of magnitude than the original ones) and the noise ξi is drawn from a Gaussian probability distribution with zero mean and variance σ2. By finding the minimizer {circumflex over (f)}j of a certain functional (to be shown in detail below), one obtains an approximation to f*j for which the average error is bounded by a constant times (k/log v)−a with 0<a≦1, i.e. the error made depends only logarithmically on v. To put it another way, the error can be made small by choosing k/log v large, but it is by no means necessary to have k/v close to 1. Accurate reconstruction is possible even if the number of projections is much smaller than v, as long as k/log v is large.
A crucial point in the derivation of the above result is the fact that the variance of the noise ξi is a constant. Even though a similar result could be obtained for non-Gaussian noises (provided certain noise properties required for the validity of the Craig-Bernstein inequality can be proved), the result does not easily carry over to the case that the variance of the noise depends on the values f*j. Yet this is precisely what happens e.g. in photon limited imaging systems where a main source of noise is the discrete quantum nature of the photons. In this case the projections yi have Poisson statistics with parameter μi=Σjφijf*j. This parameter is equal to the mean of yi but also to its variance.
In the past, it was tested whether the principle of accurate reconstructions from few projections carries over to Poisson noise in order to make accurate reconstructions possible with fewer measurements, e.g. in emission tomography. It was expected that the compressive sensing strategy for the reconstruction of sparse or compressible objects from few measurements is difficult to apply to data corrupted with Poisson noise, due to the specific properties of Poisson statistics and the fact that measurements can not usually be made randomly, as in many other cases.
Rebecca M. Willett et al. [8] have generalized results from Jarvis Haupt et al. to Poisson noise. It was proposed to reconstruct a tomographic image from detector data using a procedure of minimizing a functional {circumflex over (f)} depending on a sensing matrix A and the detector data and further depending on a penalty term, wherein the sensing matrix A is constructed on the basis of statistic Rademacher variables and the penalty term depends on the sparsity of the object. However, the result was discouraging: it was found that the upper bound on the error increases with the number of measurements, i.e. more measurements seem to make the accuracy smaller. Thus, it was assumed that compressive sensing ordinarily only works with noise types which have a fixed variance.
An image reconstruction method is proposed in [9] for reconstructing an emission tomographic image of a region of investigation within an object from detector data comprising Poisson random values measured at a plurality of different positions relative to the object. The positions are the spatial locations of measuring the detector data, e.g. the position of a detector element (pixel) of a SPECT detector device at the time of collecting data with this detector device under its current angular position with respect to the object and under its distance from the centre of rotation, or as a further example, spatial locations of detector elements sensing data at the lines of response (LOR's) of coincidence events measured with a PET detector device. According to [9], a predetermined system matrix assigning the voxels of the object to the detector data is provided. Furthermore, the tomographic image is reconstructed by minimizing a functional depending on the detector data and the system matrix and additionally including a sparse or compressive representation of the object in an orthogonal basis. The orthogonal basis is an orthogonal matrix, the columns of which are orthonormal basis vectors. The orthobasis is selected such that the object, in particular the region of investigation, can be represented in the orthobasis fulfilling the requirements of sparsity or compressibility. Contrary to the above approach of Rebecca M. Willett et al. [8] using a random sensing matrix, elements of the system matrix are not statistic values but rather selected by geometrically or physically assigning contributions of each of the voxels (object data) to each detector element of a detector device or the associated detector data, resp. The system matrix which defines the detector data as linear combinations of the original object data is determined by geometric or physical features of the imaging device, in particular by the arrangement of the object relative to the detector device and the detector device geometry.